% This script calls check to check the data and input specifications, transforms the data
% into a more easily useable form, calls estimation routine, and prints out results.
% Written by Kenneth Train on July 19, 2006, and latest edits on Sept 24, 2006.

%I moved several parts to the mxlmsl.m file.  This is just the estimation
%now.

disp('Start estimation');
disp('The negative of the log-likelihood is minimized,');
disp('which is the same as maximizing the log-likelihood.');
tic;
%parpool;
%pctRunOnAll('addpath C:\Users\Snell-PC\Box\DMNL\WTPDMIXL\AllRandom')
%parpool('AttachedFiles', {'logliknog.m','dodiscountingnog.m','llgrad2nog.m','getbetas.m'})
%pctRunOnAll('which logliknog.m')
% default algorithm is 'sqp'
%(Other available algorithms: 'active-set', 'interior-point', 'sqp-legacy', 'trust-region-reflective')
%options=optimoptions(@fmincon, 'Algorithm','sqp','MaxIterations',10000, 'OptimalityTolerance',1e-7,'StepTolerance',1e-7);
options=optimset('LargeScale','off','Display','iter','GradObj','off',...
    'MaxFunEvals',10000,'MaxIter',MAXITERS,'TolX',PARAMTOL,'TolFun',LLTOL,'DerivativeCheck','off');%, 'UseParallel', true);

A = [];
b = [];
Aeq = [];
beq = [];
nonlcon = [];
lb = [-500,-500,-500,-500,-500,-500,-500,-500, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
if discountfun=='Exponential'
    ub = [   5, 500, 500, 500, 500, 500, 500, 500, 1, 100, 100, 100, 100, 100, 100, 100, 100, 10];
elseif discountfun=="HM"
     ub = [ 5, 500, 500, 500, 500, 500, 500, 500, 2, 120, 120, 120, 120, 120, 120, 120, 120, 50];
else
      
    ub = [   5, 500, 500, 500, 500, 500, 500, 500, 2, 100, 100, 100, 100, 100, 100, 100, 100, 20];
end

if discountfun=="QuasiHyper"
    if ANA==2
        lb = [-0.0001,-50,-50,-50,-50,-50,-50,-50,-50,-50,-50,-50,-50,-50, -50,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0];
        ub = [ 1.0001, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50,  50,  1, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20];
    elseif ANA==3
        lb = [0,-10,-10,-10,-10,-10,-10,-10,-10,-10,-10,-10,-10,-10, -10,-10,-10,-10,-10,-10,-10, 0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 0];
        ub = [1, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,  10, 10, 10, 10, 10, 10, 10, 1, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 5];
    else
        lb = [-0.0001, -500,-500,-500,-500,-500,-500,-500, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
        ub = [1.0001,  500, 500, 500, 500, 500, 500, 500, 1, 100, 100, 100, 100, 100, 100, 100, 100, 10];
    end
end

printest=1;
[paramhat,fval,exitflag,output,lambda,grad,hessian]=fmincon(@logliknog,param, A, b, Aeq, beq, lb, ub, nonlcon, options);
disp(' ');
disp(' ');
%delete(gcp)
if exitflag == 1
  disp('Convergence achieved.');
elseif exitflag == 2
  disp('Convergence achieved by criterion based on change in parameters.');
  if size(PARAMTOL,1)>0
     disp(['Parameters changed less than PARAMTOL= ' num2str(PARAMTOL)]);
  else
     disp('Parameters changed less than PARAMTOL=0.000001, set by default.');
  end
  disp('You might want to check whether this is actually convergence.');

elseif exitflag == 3
  disp('Convergence achieved by criterion based on change in log-likelihood value.');
  if size(PARAMTOL,1)>0
     disp(['Log-likelihood value changed less than LLTOL= ' num2str(LLTOL)]);
  else
     disp('Log-likelihood changed less than LLTOL=0.000001, set by default.');
  end
     disp('You might want to check whether this is actually convergence.');

else
    disp('Convergence not achieved.');
    disp('The current value of the parameters and hessian');
    disp('can be accesses as variables paramhat and hessian.');
    %disp('Results are not printed because no convergence.');
    %return
end

alg_grad=grad;
alg_hessian=hessian;
if discountfun =="QuasiHyper"
    fhat = paramhat(1);
    if ANA==2
        bhat = paramhat(2:16);
        what = paramhat(17:end);
    elseif ANA==3
        bhat = paramhat(2:22);
        what = paramhat(23:end);
    else 
        bhat = paramhat(2:10);
        what = paramhat(11:end);
    end
else
    bhat = paramhat(1:9);
    what = paramhat(10:end);
end
disp(' ');
disp(['Estimation took ' num2str(toc./60) ' minutes.']);
disp(' ');
disp(['Value of the log-likelihood function at convergence/or last step: ' num2str(-fval)]);
disp(' ');
disp('FINAL Random COEFFICIENTS');
disp(' ');
disp('                Mean        Std  ');
disp('              ------------------ ');
disp('                Est         Est ');
for r=1:length(NAMES);
    fprintf('%-10s %10.4f %10.4f\n', NAMES{r,1}, [bhat(r) what(r)]);
end
disp(' ');

if NF>0
disp('FINAL Fixed COEFFICIENTS');
disp(' ');
disp('                Fixed ');
disp('              --------- ');
disp('                Est ');
for r=1:length(NAMESF);
    fprintf('%-10s %10.4f %10.4f\n', NAMESF{r,1}, [fhat(r)]);
end
disp(' ');
end
save(modelname);

disp(' ');

%Calculate standard errors of parameters
disp(' ');
disp('Taking inverse of hessian for standard errors.');
disp(' ');
ihess=inv(hessian);
stderr=sqrt(diag(ihess));
disp(stderr);
%disp(['The value of grad*inv(hessian)*grad is: ' num2str(grad'*ihess*grad)]);

%Segment parameters and account for unestimated zero mean in distribution 5
if NF>0
  fhat=paramhat(1:NF,1);
  fsd=stderr(1:NF,1);
end

if NV>0
  if sum(IDV(:,2) == 5) >0
     bhat=zeros(NV,1);
     bsd=zeros(NV,1)
     bhat(IDV(:,2) ~= 5,1)=paramhat(NF+1:NF+sum(IDV(:,2) ~= 5),1);
     bsd(IDV(:,2) ~= 5,1)=stderr(NF+1:NF+sum(IDV(:,2) ~= 5),1);
     what=paramhat(NF+sum(IDV(:,2) ~= 5)+1:end,1);
     wsd=stderr(NF+sum(IDV(:,2) ~= 5)+1:end,1);
  else
     bhat=paramhat(NF+1:NF+NV,1);
     bsd=stderr(NF+1:NF+NV,1);
     what=paramhat(NF+NV+1:NF+NV+NV,1);
     wsd=stderr(NF+NV+1:NF+NV+NV,1);
  end
end


disp('RESULTS');
disp(' ');
disp(' ')
if NF>0
disp('FIXED COEFFICIENTS');
disp(' ');
disp('                      F      ');
disp('              ------------------ ');
disp('                Est         SE ');
for r=1:length(NAMESF)
    fprintf('%-10s %10.4f %10.4f\n', NAMESF{r,1}, [fhat(r,1) fsd(r,1)]);
end
disp(' ');
end

if NV>0
disp('RANDOM COEFFICIENTS');

disp(' ');
disp('                      Mean                      SD');
disp('              ------------------   -----------------------');
disp('                 Est     SE            Est         SE');
for r=1:length(NAMES)
    fprintf('%-10s %10.4f %10.4f %10.4f %10.4f\n', NAMES{r,1}, [bhat(r,1) bsd(r,1) what(r,1) wsd(r,1)]);
end


%Create draws of coefficients from B-hat and W-hat

C=trans(bhat,what,DR);
C=reshape(C,NV,NP*NMEM);


end



disp(' ');
disp('You can access the estimated parameters as variable paramhat,');
disp('the gradient of the negative of the log-likelihood function as variable grad,');
disp('the hessian of the negative of the log-likelihood function as variable hessian,');
disp('and the inverse of the hessian as variable ihess.');
ll=-fval;
k=size(paramhat, 1);
disp('LL')
disp(ll)
disp('AIC')
disp(k*2-2*ll)
disp('BIC')
disp(k*log(NCS)-2*ll)

save(modelname);





